﻿#include <algorithm>
#include <gtest/gtest.h>
#include <iostream>
#include <vector>

#include <xy/solver/linear.h>

using namespace xy;
using namespace xy::math;
using namespace xy::solver;

#define TIMES 100

TEST(TestMatrixSolver, L)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_l::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_l solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, U)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_u::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_u solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, LU)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_lu::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_lu solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, PTLU)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_ptlu::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_ptlu solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, LLT)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_llt::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_llt solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, LDLT)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_ldlt::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_ldlt solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, Thomas)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_thomas::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_thomas solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, Blind)
{
    int count = 0;
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, A_inv, x0] = linear_blind::example(6);
        auto A2 = A;
        linear_blind solver(A2);
        solver.perform();
        auto w = solver.solve(std::move(x0));

        // 由于 w 是行向量，使用无穷范数而非 1 范数
        if (std::abs(A_inv.norm() - w.norm(0) > numeric::distance))
            count++;
    }
    // 由于 LU 分解的精度问题，有时会出现误差过大，这里只要求 80% 的测试通过
    EXPECT_LT(count, TIMES * 0.2);
}

TEST(TestMatrixSolver, Optimal)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_optimal::example(5, 6);
        auto b2 = b;

        linear_lu solver(A);
        solver.perform();
        auto x1 = solver.solve(std::move(b2));

        linear_optimal optimal(A, x1);
        optimal.perform();
        auto x2 = optimal.solve(std::move(b));

        // 迭代后 x2 与 x1 相比，误差应该小于 x1 误差的 10 倍
        A.transpose();
        auto r1 = (x - x1).norm();
        auto r2 = (x - x2).norm();
        EXPECT_LT(r2, r1 * 10);
    }
}

TEST(TestMatrixSolver, Jacobi)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_jacobi::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_jacobi solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, GaussSeidel)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_gauss_seidel::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_gauss_seidel solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, SOR)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_sor::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_sor solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, QR)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto [A, x, b] = linear_qr::example(5, 6);
        auto A2 = A;
        auto b2 = b;
        linear_qr solver(A2);
        solver.perform();
        auto x2 = solver.solve(std::move(b2));
        EXPECT_LT((x - x2).norm(), numeric::distance);
    }
}

TEST(TestMatrixSolver, EigenVector)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto A = matxd::symmetric_positive_definite(5);
        linear_eigen solver(A);
        solver.perform();
        auto lambda = eigenvalue(A);
        auto v = solver.solve(lambda).transposed();
        // std::cout << solver.move_matrix() << std::endl;
        EXPECT_LT((A * v - v * lambda).norm(), 1e-1);
    }
}

TEST(TestMatrixSolver, PassJacobi)
{
    for (int i = 0; i < TIMES; i++)
    {
        auto A = matxd::symmetric_positive_definite(5);
        linear_sym_eigen solver(A);
        solver.perform();
        auto Q = solver.solve(10);
        // std::cout << solver.move_matrix() << std::endl;
        EXPECT_LT((Q.transposed() * A * Q - solver.move_matrix()).norm(), numeric::distance);
    }
}